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Abstract. We analyze stress transmission in wet granular media in the pendular state by means of three- 
dimensional molecular dynamics simulations. We show that the tensile action of capillary bonds induces 
a self-stressed particle network organized in two percolating "phases" of positive and negative particle 
pressures. Various statistical descriptors of the microstructure and bond force network are used to charac- 
terize this partition. Two basic properties emerge: 1) The highest particle pressure is located in the bulk 
of each phase; 2) The lowest pressure level occurs at the interface between the two phases, involving also 
the largest connectivity of the particles via tensile and compressive bonds. When a confining pressure is 
applied, the number of tensile bonds falls off and the negative phase breaks into aggregates and isolated 
sites. 
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1 Introduction 

The particle-scale origins of the strength and flow prop- 
erties of dry granular materials has been a subject of in- 
tensive research since several decades. It is now generally 
admitted that the scale-up of particle interactions to the 
macroscopic scale is more subtle than initially expected 
because of mediation by a disordered microstructure with 
a rich statistical behavior |H2j . Considerable work has 
thus been devoted to characterize the microstructure and 
its manifestations in the form of highly inhomogeneous 
distributions of interparticle forces and fluctuating parti- 
cle velocities [314 5 6 7 8 91 10111112113] . The bimodal char- 
acter of the force network [II], exponential probability 
functions of strong forces (force chains) 15 16 8117] . and 
clustering of dissipative contacts are recent examples of 
this non trivial phenomenology [18j . Insightful analogies 
have also emerged with other fields such as jamming tran- 
sition in colloidal matter |19[20|21|22j . slow relaxation in 
glassy materials [23], and fluid turbulence [12] . 

Most of our present knowledge on the subject excludes, 
however, cohesive bonding between particles. Although 
one expects strong similarities due to the common granu- 
lar microstructure, the presence of cohesion brings about 
new mechanisms that tend to transform the nature of 
the problem. At the macroscopic scale, the shear strength 
needs to be described in terms of the Coulomb cohesion in 
addition to the internal angle friction [24]. On the other 
hand, the interplay between cohesive bonds, friction and 
rotation frustration of the particles leads to novel features 
such as particle aggregation that control static and dy- 
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namic properties of the material ,25^26]. A well-known 
example is the wet sand where small amounts of water 
affect significantly the bulk behavior [27128129] . The phe- 
nomena arising from cohesion are of particular interest to 
the processing (compaction, granulation. . . ) of fine pow- 
ders [30131] . It seems thus that a systematic investigation 
of the microstructure in cohesive granular media should 
open new scopes for modeling granular materials. 

In this paper, we analyze the force network in wet 
granular assemblies of spherical particles. We are inter- 
ested in a basic question: how cohesive grains keep to- 
gether to form a self-sustained structure in the absence 
of confinement (no container and no confining stresses)? 
The packing can reach an equilibrium state as a result of 
attraction forces and elastic repulsion between particles 
without or with self-stressed structures. While the parti- 
cles are balanced in both cases, the attraction force at each 
contact is exactly balanced by an elastic repulsion force in 
the first case. In contrast, in the second case all contacts 
are not in their equilibrium state due to steric hinderance 
between particles. Hence, a network of tensile and com- 
pressive forces is formed inside the packing. These self- 
equilibrated forces can be induced through various loading 
histories such as consolidation [32 33 or differential par- 
ticle swelling 34 . In wet granular media in the pendular 
state, the self-stresses appear naturally due to the tensile 
action of capillary bonds bridging the gaps between neigh- 
boring particles within a debonding distance. We focus in 
this paper on the structure of these self-stresses induced 
by capillary bonds. 

We use 3D molecular dynamics method in which capil- 
lary attraction between spherical particles is implemented 
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as a force law expressing the capillary force as a function 
of the distance, water volume, and particle diameters. The 
total water volume is distributed homogeneously between 
particles. The packing is analyzed at equilibrium under 
zero confining pressure. The main goal of our analysis 
is to characterize the organization of particle pressures 
which take positive or negative values according to their 
positions in the network of self-equilibrated bond forces. 
We will see that this organization involves a genuin parti- 
tion of the particles in two phases of negative and positive 
pressures. In the following, we first describe the simulation 
method and our model of capillary cohesion. In Sections [3] 
and|H we study in detail the probability density functions 
of the forces and the connectivity of particles via tensile 
and compressive bonds. Then, in Sections [5] and [6] we in- 
troduce particle stresses and we analyze their distribution 
and correlation with the connectivity of the particles. Sec- 
tion [7] is devoted to the influence of external pressure. We 
conclude with a discussion about the main findings of this 
work. 



2 Numerical method 

We used the molecular dynamics (MD) method with a 
velocity Verlet integration scheme [35 36]. The force laws 
involve normal repulsion, capillary cohesion, Coulomb fric- 
tion, and normal damping. The normal force has three 
different sources, 



fn — fn + fn + fn' 



(1) 



The first term is the repulsive contact force depending 
linearly on the normal distance 5 n between the particles 
(FigfUa)): 
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where k n is the normal stiffness. The second term repre- 
sents a viscous damping force depending on the normal 
velocity S n : 



J n 
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where m = mirrij / (rrii + rrij) is the reduced mass of the 
particles i and j, a n is a damping rate varying in the range 
[0, 1[ and that accounts for the rate of normal dissipation 
or the restitution coefficient between particles [57] . 

The last term in Eq. [T]is the capillary force depending 
on the liquid bond parameters, namely the gap 5 n , the 
liquid bond volume VJ,, the liquid surface tension j s , and 
the particle-liquid-gas contact angle 9. The capillary force 
can be obtained by integrating the Laplace- Young equa- 
tion [38 39 40 41 . However, for molecular dynamics simu- 
lations, we need an explicit expression of f° as a function 
of the liquid bond parameters. We propose the following 
simple form: 
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Fig. 1. (a) Geometry of a capillary bridge; (b) Capillary force 
as a function of the gap 8 n between two particles for different 
values of the liquid volume Vb and size ratio r according to the 
model proposed in this paper (solid lines), and from direct 
integration of the Laplace- Young equation (open circles); (c) 
Scaled plot of the capillary force as a function of gap from the 
direct data shown in (b). 



where R = */ RiRj is the geometrical mean of particle 
radii, and A is a length scale to be discussed below. The 
prefactor k is given by [42 43 44J 



k = 2tt^/ s cos 9, 
and 5™ aa; is the debonding distance given by 



(5) 



(6) 
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The capillary bridge is stable as long as S n < 5™ ax . In the 
simulations, the bridge is removed as soon as the debond- 
ing distance is reached, and the liquid is redistributed 
among the contacts belonging to the same particle in pro- 
portion to grain sizes [35] • At contact (S n — 0), kR corre- 
sponds to the largest attraction force between two parti- 
cles. In the simulations presented in this paper, we assume 
that particles are perfectly wettable, i.e. 8 = 0. This is a 
good approximation for water and glass beads. 

The length A governs the exponential falloff of the cap- 
illary force in Eq.[U It should depend on the liquid volume 
Vb, amean radius R' , and the ratio r = max{i?j/i?j ; Rj/Ri} 
A reasonable choice is 

\ = ch(r)(^y\ (7) 

where c is a constant prefactor, h is a function depend- 
ing on the ratio r, and R' is the harmonic mean (R' = 
2RiRj/(Ri+Rj)). When introduced in Eqs.[7]and[31 this 
form yields a nice fit for the capillary force obtained from 
direct integration of the Laplace- Young equation by sim- 
ply setting h(r) = r^ 1 / 2 and c ~ 0.9. Fig. QJb) shows 
the plots of Eq. [4] for different values of the liquid volume 
Vb and size ratio r together with the corresponding data 
from direct integration. The fit is excellent for 8 n = (at 
contact) and for small gaps. 

Figure [2(c) shows the same plots of the direct data as 
in Fig. [2(b) but in dimensionless units where the forces 
are normalized by kR and the lengths by A. We see that 
the data collapse indeed on the same plot, indicating again 
that the force kR and the expression of A in Eq. [7] char- 
acterize correctly the behavior of the capillary bridge for 
8 = 0. From the same data, we also checked that the geo- 
metric mean R = ^jRiRj introduced in Eq. [4] provides a 
better fit than the harmonic mean QR^Rj j [R,% -t- Rj) pro- 
posed by Derjaguin for polydisperse particles in the limit 
of small gaps [46] . 

For the friction force ft, we used a viscous-regularized 
Coulomb law [37147148] 

ft = -min h\\S t \\, M (/ n -/=)}**-, (8) 

1 } INI 

where 7* is a tangential viscosity parameter, (i is the co- 
efficient of friction, and dt is the sliding velocity. In relax- 
ation to equilibrium, S t declines but never vanishes due 
to residual kinetic energy. The equilibrium state is prac- 
tically reached as soon as we have 7t||5t|| < /i(/„ — 
at all contacts. This means that the friction force is inside 
the Coulomb cone everywhere in the system. 

The simulations were performed with a packing com- 
posed of N = 8000 spheres with diameters d — 1, 1.5, and 
2 mm, in equal numbers. The system was subjected to an 
isotropic pressure p m via six rigid walls and no gravity in 
order to obtain an as homogeneous state as possible. For 
the same reason, the friction with the walls was set to zero 
although wall effects can not be fully removed. The param- 
eters used in the simulations are displayed in Table [TJ The 
choice of the water volume has no influence on the value of 



Table 1. Parameters used in the simulations. 



Parameter 


Symbol 


Value 


Unit 


Normal stiffness 


fen 


1000 


N/m 


Damping rate 


Oin 


0.8 




Tangential viscosity 


It 


1 


Ns/m 


Capillary force prefactor 


K 


0.4 


N/m 


Gravimetric water content 




0.007 




Liquid density 




1000 


kg-.m" 3 


Particle density 




2700 


kg.m 


Friction coefficient 




0.4 




Time step 




10" 6 


s 



the largest capillary force in the pendular state. We also 
note that the liquid bonds are homogeneously distributed 
into all gaps within the debonding distance [39]. With a 
gravimetric water content of 0.007, the coordination num- 
ber is ~ 6. Experiments show that the distribution of liq- 
uid bonds depends on the preparation protocol involving 
the water volume, mixing procedure, and waiting times 
[501391 . 



3 Force distributions 

We start out by considering force distributions first in a 
system subjected to a negligibly small average compressive 
stress p m ~ Pa. Fig. [3J displays the probability density 
function (pdf) of the normal forces both in tensile (nega- 
tive) and compressive (positive) ranges. We observe a dis- 
tinct peak centered on f n = and two nearly symmetrical 
parts decaying exponentially from the center: 

P oce-^l/A, (9) 

with a ~ 4 within statistical precision for both negative 
and positive forces, and /q = nR max , where R m ax is the 
largest particle radius. It has been observed that in dry 
granular media the distribution is not purely exponen- 
tial in the whole range of bond forces [16 8 I17I5T] . Below 
the average force, the distribution tends to be uniform 
or a decreasing power law with a weak exponent [55] . In 
the present case of wet cohesive materials, the exponen- 
tial behavior extends to the center of the distribution. It 
is important to remark that this peak does not represent 
non-transmitting contacts. It rather corresponds to con- 
tacts where capillary attraction is balanced by elastic re- 
pulsion, i.e. k n 8 n + fo = 0. We will see below that this 
peak persists under the action of a compressive confin- 
ing stress, suggesting that its presence reflects a feature of 
force transmission in wet granular materials. 

The tensile range is cut off at /„ = — fa corresponding 
to the largest capillary force. Although the confining stress 
is zero, positive forces as large as 2/ can be found in 
the system. In contrast to dry granular materials, the pdf 
shows a peak at f n = which is the average force in 
the present case. In fact, in an unconfined assembly of 
dry rigid particles, no self-stresses appear and the forces 
vanish at all contacts. In our wet material, the presence of 
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Fig. 2. Probability density function of normal forces normal- 
ized by the largest capillary force fo for zero confining pressure. 




Fig. 3. (Color online) A map of tensile (green) and compressive 
(red) forces in a thin layer cut out in the packing. Line thickness 
is proportional to the magnitude of the force. 

liquid bonds induces tensile and compressive self-stresses 
although the average force is zero. 

Figure [3] shows the force network in a narrow slice 
nearly three particle diameters thick. The tensile and com- 
pressive forces are represented by segments of different 
colors joining particle centers. The line thickness is pro- 
portional to the force. It is remarkable that tensile and 
compressive force chains can be observed although the 
slice is quite narrow. The bond coordination number z 
(average number of bonds per particle) is ~ 6.1 including 
nearly 2.97 compressive bonds and 3.13 tensile bonds. 

4 Connectivity of the bond network 

We analyze the connectivity of the particles via capillary 
bonds by considering the fraction p{k + , k~) of particles 



k~ 




Fig. 4. Grey level map of the connectivity function p(k + , k ). 

with exactly k + compressive bonds and k~ tensile bonds. 
This function is displayed in Fig. [4] as grey level map in 
the parameter space (k + , k~) for our system. The map 
is symmetric with respect to the line k + — k~ , reflecting 
thus the symmetrical roles of compressive and tensile net- 
works in the absence of confining stresses. A peak value 
of p occurs at fc + = k~ = 2. Obviously, the condition 
of particle equilibrium cannot be fulfilled with k + < 1 
and k~ < 1 and the corresponding levels are zero on the 
map. The particles with k + = 2 and k~ = define chains 
of compressive bonds whereas particles with k + — and 
k~ = 2 correspond to chains of tensile bonds. But such 
"pure" chains are rare. At larger values of k + and k~ the 
fractions decline basically due to geometrical hinderance 
between particles. 

An interesting feature of the connectivity map (Fig. 2]) 
is the existence of a population of particles involving no 
tensile bonds (the row k~ = 0) as well as a population 
of particles involving no compressive bonds (the column 
k + = 0). Reduced connectivity functions p + and p~ can 
be defined by summing up the function p(k + , fc~) along 
columns and rows, respectively: 

P + (k + ) = Y / P(k + ,k~), 

fe- 

p-(k-) = Y.P( k+ > k ~)- ( 10 ) 

k+ 

The plots of these functions are shown in Fig. [5j They 
nearly coincide as expected from the symmetry observed 
in Fig. |U We have p + (0) = p~(0) — 0.08, corresponding 
respectively to particles in purely extensional or compres- 
sional local environments. A maximum occurs at k + = 
k~ = 2 or 3. 



5 Particle stresses 

Until now, we focussed on forces and their distributions 
with regard to the tensile or compressive nature of the 
bonds. For the description of stress transmission in our 
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Fig. 5. Reduced connectivity functions at zero confining pres- 
sure for tensile and compressive bonds. 
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Fig. 6. Probability density function of particle pressures nor- 
malized by reference pressure po (see text) in the unconfincd 
packing. 



system we need, however, to characterize the inhomogencitics 
at the scale of particles representing the smallest entities 
for which the equations of motion are resolved in MD 
simulations. At this scale, a Cauchy stress in the sense 
of continuous media cannot be defined. But, it can be 
shown that the so-called internal moment tensor M has 
the same properties as the Cauchy stress tensor er, and 
its definition extends mathematically to discrete media at 
all scales down to the particle scale |53|54j . For a particle 
i subjected to forces P J from its contact neighbors j at 
the contact points r u , the internal moment tensor Mj is 
given by [53] 

A/, r'. (11) 

where ® designs a tensorial product. The internal moment 
tensor is additive and independent of the origin of the 
coordinate frame. For a collection of particles in a control 
volume V, the total internal moment M is simply the sum 
of the particle moments: 



following form: 



(14) 



(12) 



This tensor has the dimension of a moment and it becomes 
homogeneous to a stress when divided by the control vol- 
ume: 

' T = T7E M - ( 13 ) 

i£V 

At large scales (containing a sufficiently large number of 
particles), the volume is well-defined and the stress tensor 
is equivalent to the internal moment tensor divided by 
this volume. However, at the particle scale, while Mi is 
defined in a unique way, the choice of the volume V is a 
matter of convenience. It is reasonable to take into account 
the free volume Vi of each particle i, as the sum of the 
volume of the particle and part of the surrounding pore 
space. On average, we have Vi = (l/6)ndf /v, where di is 
the particle diameter and v denotes the packing fraction. 
With this choice, the particle stress tensor <?i takes the 



Remark that when summing this form over many particles 
in a control volume as in Eq. I13( each contact ij enters 
the summation two times, belonging once to particle i and 
once to particle j. Since f u = — f J \ the contribution of the 
contact ij to stress is P J ®(r y — r Ji ), a writing that involves 
the center-to-center vector instead of position vectors of 
the contact points if the origin of coordinates for each 
particle is chosen coincides with the center of the particle. 
However, we consider here only the particle stress, and at 
this scale, according to Eq. [TU only the position vectors 
of contact points are involved. 

Here, we explore the particle pressures (average par- 
ticle stresses) pi = (l/3)tr<Xj. Each particle can take on 
positive or negative pressures according to the nature of 
the forces exerted by contact neighbors. The pdf of par- 
ticle pressures is displayed in Fig. The pressures have 
been normalized by a reference pressure po = fo/(d} 2 . 
The distribution is symmetric around and peaked on zero 
pressure, and each part is well fit by an exponential form. 
Obviously, the exponential shape of particle pressure dis- 
tributions reflects statistically that of bond forces. In dry 
granular media, since the normal forces are all of the same 
sign (compressive) and particle pressure results from the 
summation of individual bond forces on a particle, the 
probability is expected to vanish as the pressure goes to 
zero. In contrast, Fig. [6]shows that the exponential shape 
of particle pressure distribution extends to the center of 
the pdf. This can be related to the fact that all normal 
forces are not of the same sign. On the other hand, the 
zero pressure corresponds to a state where a particle is 
balanced under the combined action of tensile and com- 
pressive forces. Since such particle states are not marginal 
here, they might reflect a particular organization of the 
stresses in the wet packing (see below). 
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Fig. 7. (Color online) The map of partial pressures p(k^ 

6 Partition of particle pressures 



k~). 



An interesting observation about the connectivity map 
(Fig. [4]) was the presence of particles with only tensile 
or only compressive bonds. In terms of particle pressures, 
these populations carry negative or positive pressures, re- 
spectively. However, this information is not rich enough 
as it does not specify the pressure levels in these popula- 
tions. The question is how the particle pressure is locally 
correlated with the particle connectivity. For particle i, 
the connectivity is specified by the number kf~ of com- 
pressive bonds and the number fer of tensile bonds. Let 
us now consider the set S(k + , k~) of particles i such that 
kf = k + and fc~ = k~. The partial pressure carried by 
this set is the sum of particle pressures in this set divided 
by the total number of particles: 



p(k + , k~) 



1 E 



Pi- 



(15) 



ieS(k+,k-) 



This is obviously an additive quantity so that the average 
stress p m = J2(k+ k-)P(k + i k~). The function p(k + , k~) 
provides detailed information about the way particle pres- 
sures are distributed with respect to the bond network. In 
other words, this function describes the relationship be- 
tween the pressure sustained by a particle and its first 
neighbors with which it is bonded. 

Figure [7] shows the map of partial pressures p(k + , fc~) 
in the parameter space (k + , k~). Interestingly, we observe 
a bipolar structure of partial pressures which is antisym- 
metric with respect to the line k + = fc~ within statis- 
tical precision. The pressures are positive in the range 
k + > fc~ and negative in the range k + < k~ . Hence, 
the line k + = k~ defines the transition zone between 
the two parts with p(k + , k~ = k + ) ~ 0. This line cor- 
responds to the largest population of particles according 
to the connectivity map (Fig. [4j. The pressure extrema 
are located at (k + = 4, k~ = 0) for positive pressures and 
at (k + = 0, k~ = 4) for negative pressures. 

The bipolar structure of the pressure map suggests 
that the particles of positive and negative pressures define 




Fig. 8. A representation of the packing with particles of neg- 
ative (white) and positive (black) pressures. 




Fig. 9. Average numbers of tensile (z~) and compressive (z + ) 
bonds per particle as well as the partial coordination number 
(z — z~ + z + ) as a function of the particle pressure in the 
unconfined packing. 



two separate phases throughout the system. In this pic- 
ture, the population of particles along the line k + = k~ 
corresponds to the particles at the interface between the 
two phases. This interpretation is backed by Fig. [5] dis- 
playing the packing where the two phases are represented 
in black and white. The particles of either positive or neg- 
ative pressure percolate throughout the system. The two 
phases are intimately intermingled with a large interface 
between them. The particles at the interface belong to the 
line k + = fc~ corresponding to a fraction 0.125 of parti- 
cles. The morphology of the two phases is approximately 
filamentary with variable thickness. 

The correlation between the bond connectivity and 
particle pressures can alternatively be determined by con- 
sidering the average numbers of tensile and compressive 
bonds per particle, denoted z~ and function of the 

particle pressure p. In order to evaluate these functions, 
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we consider the set S(p) of particles i with a pressure pi in 
the range [p, p + Ap], where Ap is the pressure increment. 
The partial coordination numbers z~ and z + arc dehncd 
by 



(16) 



where iV(p) is the number of particles in the set. These 
functions arc plotted in Fig.[S]in the case p m — 0. The plot 
of z~(p) is an approximate mirror image of z + {p) with 
respect to p = 0. Three zones can be discerned. For p < 
—po, we have z~ ~ 5 and z + ~ 0.5. This zone represents 
mainly the particles belonging to the bulk of the negative 
phase. For p > p , we have z~ ~ 0.5 and z + ~ 4. The 
particles in this zone belong to the bulk of the positive 
phase. In the range — po < p < pa, z + increases and z~ 
declines. The intersection occurs at p — where z~ = 
z + ~ 3. This zone corresponds to the particles located at 
the interface between negative and positive phases. 

The above findings underline that stress transmission 
in wet granular media is non trivial. In particular, they 
support the partition of the packing into two well-defined 
phases both in terms of particle connectivities and in terms 
of particle pressures. The peculiarity of this partition is 
that the extrema of particle pressures are located in the 
bulk of each phase (Fig. |4]) whereas the maximum of con- 
nectivity between particles via tensile and compressive 
bonds resides at their interface (Fig. EJ. Along this in- 
terface, not only the bond forces are balanced on each 
particle, as everywhere in the packing, but the tensile and 
compressive bonds contribute equally in such a manner 
that the particle pressures are extremely low. 



7 Influence of confining pressure 

In the absence of a confining stress, the capillary bonds are 
at the origin of self-stresses or self-equilibrated forces that 
we analyzed in preceding sections. Hence, the observed 
symmetry between tensile and compressive bonds is a con- 
sequence of static equilibrium at zero confining stress. The 
question remains whether the partition of particle pres- 
sures, as depicted above, still holds when a wet packing 
is subjected to (compressive) confining stress. In practice, 
however, we cannot isolate self-equilibrated forces for each 
particle from those induced by the external stress (the ac- 
tual force network being the sum of the two networks). 
This is because the external pressure drives the pack- 
ing to a new equilibrium state with modified microstruc- 
ture. Hence, the self-stresses for p m — 0, i.e. before rear- 
rangements, do not correspond to the rearranged state for 
p m 7^ 0. While it can be conjectured that the self-stresses 
in the presence of a confining pressure will display the 
same "bipolar" behavior as for p m = 0, we consider here 
the force network and particle pressures without distinc- 
tion between induced and self-equilibrated forces. 
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Fig. 10. Probability density function of normal forces normal- 
ized by the largest capillary force fo in the confined packing. 



We applied an isotropic stress p m = 100 Pa to the 
wet packing prepared with p m — 0. The packing was then 
allowed to relax to equilibrium under the action of the ap- 
plied pressure. This level of confinement is high compared 
to the reference pressure po (p m /po — 0.5), yet not too 
high to mask fully the manifestations of capillary cohe- 
sion. The same packing was also compressed isotropically 
forp m = 100 Pa without capillary cohesion (dry packing). 

The pdf of normal forces is shown in Fig. [10] for both 
packings. The symmetry of the distribution around /„ = 
is now broken (compare to Fig. [5]). The distributions 
are roughly exponential for both tensile and compressive 
forces but the exponents are different. We also note that 
the exponent for compressive forces is nearly the same 
as for the dry sample. On the other hand, although the 
strictly zero forces have been removed from the statistics, 
a distinct peak centered on zero force is present for the wet 
sample and absent from the dry sample. If the occurrence 
of this peak in the unconfined case (p m = 0) is attributed 
to the interfacial zone, its persistence in the confined case 
suggests that a negative pressure phase continues to co- 
exist with the positive pressure phase (which is now the 
dominant phase assisted with the confining stress). This 
point will be analyzed below in relation to the distribution 
of particle pressures. 

The coordination numbers for compressive and ten- 
sile bonds are 4.85 and 0.85, respectively. This shows that 
the application of a confining pressure has transformed 
a fraction of tensile bonds into compressive bonds. The 
plots of the connectivity functions p + (k + ) and p~(k~) 
are displayed in Fig. 1111 Each function is normalized to 
unity as in Fig. [5] for the unconfined packing. The func- 
tion p~(k~) decreases with k~ from a peak at fc~ = 0, 
showing that nearly half of the particles have no tensile 
contact at all. The function p + {k + ) has a peak at fc + = 4 
and only a small fraction ~ 0.05 of particles have no com- 
pressive bond (p + (k + = 0) ~ 0.05). 

Force patterns with one or more tensile bonds corre- 
spond to various local configurations where equilibrium of 
the particles cannot be ensured only by compressive con- 
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Fig. 11. Reduced connectivity functions in the confined pack- 
ing for tensile and compressive bonds. 




Fig. 12. (Color online) Two examples of local patterns of ten- 
sile (green) and compressive (red) forces surrounding particles 
with negative (green) and positive (red) pressures. 



tacts. Typically, a tensile bond between two particles is 
induced by transverse particles that are forced into the 
space between the two particles. One example is shown in 
Fig. HHa). For this reason, when a packing is subjected 
to a stress deviator, most tensile bonds occur along the 
direction of extension 29J . The upshot of tensile bonds 
when they are located in a purely compressive environ- 
ment is thus to reinforce the shear strength. We also ob- 
serve many particles with 2, 3 or 4 tensile bonds. One 
example is shown in Fig. I12f b) where the tensile actions 
of several particles creates a "cage" of compressive bonds 
between their neighboring particles. This kind of patterns 
may be locally self-equilibrated and thus form aggregates 
that could move as a rigid body when the packing deforms. 

The pdf of particle pressures is shown for dry and wet 
samples in Fig. [T3l Large positive pressures decay expo- 
nentially in both cases. In the dry case, a local maximum 
is observed at p ~ l-5po as a signature of the confining 
stress. In this range of pressures, we observe a plateau for 
the wet sample. In the range of negative pressures, the 
distribution is no more exponential. This means that the 
organization of tensile forces does not fulfill the conditions 
that underly exponential distribution of strong forces in 
granular media |15j . In particular, the network of tensile 
forces is no more percolating throughout the packing as in 
the unconfined case. A map of positive and negative par- 
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Fig. 13. Probability density function of particle pressures nor- 
malized by a reference pressure po (see text) in the confined 
packing. 




Fig. 14. A representation of a thin layer inside the confined 
packing with negative (white) and positive (black) particle 
pressures. 



ticle pressures in a thin layer inside the packing is shown 
in Fig. 1141 The positive pressures are dominant and nega- 
tive pressure particles are mostly isolated or appear in the 
form of small clusters. Although the negative pressures do 
not define a bulk phase any more, the peak centered on 
zero pressure in Fig. [T3] can still be considered as a remi- 
niscence of the interface between the two phases. 

This last point appears more clearly in the plots of 
partial coordination numbers z~ and z + as a function of 
the particle pressure p in Fig. [151 We again discern three 
zones as in Fig. |H] for the unconfined packing. The peak 
z + ~ 6 appearing around p ~ 2.5 po is the effect of the 
confining pressure. At the same time, the level of z~ in 
the range of negative pressures is reduced to ~ 3.5 (from 
~ 5 in the unconfined packing). The intersection occurs 
at p = with z~ — z + ~ 2. 
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Fig. 15. Average numbers of tensile {z~) and compressive (z + ) 
bonds per particle as well as the partial coordination number 
(z — z~ + z + ) as a function of the particle pressure in the 
confined packing. 




Fig. 16. A representation of the packing with particle pres- 
sures above (white) and below (black) the m6£Lii pressure pm • 

We see that many features of stress transmission in the 
unconfined packing persist when a confining stress is ap- 
plied, fn particular, in both large class of particles 
of weak pressure (close to zero of either sign) is present. 
This class was interpreted in the unconfined case as be- 
longing to the interface between two percolating phases. 
Obviously, this interface is ill-defined for p m = 100 Pa 
where the negative phase appears in the form of either 
isolated particles or very small aggregates. Nevertheless, 
our results clearly indicate that tensile bonds and negative 
pressures play the same role with respect to the equilib- 
rium properties of the particles wherever they are present. 

Obviously, as stated before, the molecular dynamics 
method cannot be used as such for the analysis of self- 
stresses in the presence of an external pressure. This is 
because in molecular dynamics, the resolution of govern- 
ing equations proceeds from the knowledge of the posi- 



tions of only the first neighbors of each particle and not 
from an explicit construction of the system of equations 
as in contact dynamics [55156] An approach for the analy- 
sis of self-stresses was presented in Radja'i et al. [57] using 
"singular value decomposition" in the framework of the 
contact dynamics method. However, a rough estimation 
of the self-stresses can be obtained by simply subtracting 
the mean pressure p m from the particle pressures. Fig. 1161 
displays a snapshot of the particle pressures where the 
pressures below and above p m are distinguished. The ap- 
parent clustering of particle pressures is quite comparable 
to that observed in Fig [8] This indicates that it is very 
likely that if the self-stresses were isolated from those in- 
duced by the external pressure, they would display the 
same nearly symmetric "bipolar" structure as that ob- 
served at zero confining pressure. 



8 Conclusions 

We analyzed the statistical properties of the network of 
self-equilibrated forces in a wet granular material by means 
of 3D molecular dynamics simulations. Various descrip- 
tors of the microstructure and bond force network were 
shown to carry the signature of an ingenious organization 
of particle pressures in two distinct clusters of respectively 
positive and negative pressures, each percolating through- 
out the packing. This partition is not meant as a formal 
distinction between negative and positive pressures but 
rather related to the way the two populations share the 
space and connect to the bond network. This "phase sep- 
aration" is characterized by two interesting properties. 
First, the highest pressures occur in the heart of each 
phase, whereas the lowest pressure levels constitute the 
interface between the two phases. Secondly, this interface 
bears the largest coordination numbers via tensile and 
compressive bonds. In the presence of confining stresses, 
the same phenomenology can be expected for self-stresses 
although these can not be directly accessed from the force 
data. 

It is important to remark that the homogeneity of self- 
stresses in our simulations results from the homogeneous 
distribution of capillary bonds. Obviously, the self-stresses 
can be more or less localized in different portions of the 
material or involve gradients if the liquid bonds are dis- 
tributed in a nonuniform manner in the bulk. In the same 
way, although the boundary conditions are isotropic, the 
self-stresses may be locked in an anisotropic state as a re- 
sult of friction and geometrical hinderance effects. In par- 
ticular, if the capillary bonds are placed only in the gaps 
between particle pairs with privileged orientation, the self- 
stresses might organize into an anisotropic scheme. Such 
anisotropic distributions of liquid bonds may also appear 
as a result of handling the material. The choice of a ho- 
mogeneous distribution of liquid bonds in our simulations 
was motivated by the requirement of representative statis- 
tics for the analysis of the system. However, it would be 
interesting to study the influence of the liquid distribution 
on the patterns of self-stresses. 
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The partition of self-stresses implies that the overall 
equilibrium of the packing is ensured by mesoscopic struc- 
tures involving length scales larger than the particle size. 
These length scales are likely to control the size of aggre- 
gates during flow or other packing properties of cohesive 
granular materials. On the other hand, the effect of self- 
stresses on the tensile strength or Coulomb cohesion of 
wet granular materials is of interest to wet processing of 
grains in chemical engineering and merits to be studied 
along these lines, fn the same way, the influence of solid 
fraction is an important aspect with evident application 
to compaction and consolidation of cohesive packs. 

This work was accomplished within the "granular solids" group 
of the LMGC. We thank F. Soulie for his help with the valida- 
tion of the capillary law used in our simulations. The data were 
treated by means of the 3D software mgpost (www.lmgc.univ- 
montp2.fr/~richefeu). 
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